/*
 * Interpolation3DBase has the following inherited classes:
 * - PolInt3D for General Polynomial Interpolation for 3 dim
 * - LinInt3D for Efficient Linear Interpolation for 3 dim
 * - SplInt3D for Spline Interpolation for 3 dim
 * 
 * Author: Yuejia Zhang
 * Date: 2023/2/21
 * 
 * */

#ifndef INTERPOLATION_3D
#define INTERPOLATION_3D

#include <cmath>
#include <iostream>

#include "Interpolation2D.h"

template <class InternalType = double, class ExternalType = float>
class Interpolation3DBase {
 protected:
  int m_, n_, l_;  // number of intervals, one less than number of knots
  ExternalType *xa_, *ya_, *za_;
  ExternalType ***va_;
  InternalType *coefficients;

 public:
  /* The construct function computes and stores the coefficients
   * (or second derivatives). xa, ya, za are arrays of x, y, z,
   * va is 3d-array of v. m is the dimension of xa,
   * n is the dimension of ya, l is the dimension of za. */
  virtual void construct(ExternalType xa[], ExternalType ya[],
                         ExternalType za[], ExternalType ***va) = 0;

  /* The evaluate function requires to find the place in the table,
   * at worst O(log N). */
  virtual ExternalType evaluate(ExternalType x, ExternalType y,
                                ExternalType z) = 0;

  /* subfunction of evaluate, using binary search to locate the interval
   * (xa[i], xa[i+1]) where x is in. returns the left index i.
   * xa must be a sorted array in increasing order!
   * also handles the boundary condition, see the first two lines. */
  int binarySearch(ExternalType *xa, int n, ExternalType x) {
    if (x < xa[0]) return -1;
    if (x > xa[n - 1]) return n;
    int ib = 0, ie = n - 1;
    int im = (ib + ie) / 2;
    while (ib < im) {
      if (x < xa[im])
        ie = im;
      else
        ib = im;
      im = (ib + ie) / 2;
    }
    return im;
  }

  /* subfunction of evaluate, use linear interpolation when either
   * x < x_min or y < y_min or z < z_min
   * the linear function is induced from grid
   * (xa[ib], xa[ib+1]) x (ya[jb], ya[jb+1]) x (za[kb], za[kb+1]). */
  ExternalType evaluateLinear(ExternalType x, ExternalType y, ExternalType z,
                              int ib, int jb, int kb) {
    ExternalType *xa = this->xa_, *ya = this->ya_, *za = this->za_;
    ExternalType ***va = this->va_;
    InternalType dx = xa[ib + 1] - xa[ib], dy = ya[jb + 1] - ya[jb],
                 dz = za[kb + 1] - za[kb], xdx = x - xa[ib], ydy = y - ya[jb],
                 zdz = z - za[kb];
    InternalType c000, c001, c010, c011, c100, c101, c110, c111, dt;
    c000 = va[ib][jb][kb];
    c100 = (va[ib + 1][jb][kb] - va[ib][jb][kb]) / dx;
    c010 = (va[ib][jb + 1][kb] - va[ib][jb][kb]) / dy;
    c001 = (va[ib][jb][kb + 1] - va[ib][jb][kb]) / dz;
    c110 = ((va[ib + 1][jb + 1][kb] - va[ib][jb + 1][kb]) / dx - c100) / dy;
    c011 = ((va[ib][jb + 1][kb + 1] - va[ib][jb][kb + 1]) / dy - c010) / dz;
    dt = (va[ib + 1][jb][kb + 1] - va[ib][jb][kb + 1]) / dx;
    c101 = (dt - c100) / dz;
    c111 = (((va[ib + 1][jb + 1][kb + 1] - va[ib][jb + 1][kb + 1]) / dx - dt) /
                dy -
            c110) /
           dz;
    return c000 + c100 * xdx + (c010 + c110 * xdx) * ydy +
           (c001 + c101 * xdx + (c011 + c111 * xdx) * ydy) * zdz;
  }

  /* help function for solely debug purpose. */
  virtual void printCoefficients() = 0;
};

/* LinInt3D: Efficient Linear Interpolation for 3 dim
 * */
template <class InternalType = double, class ExternalType = float>
class LinInt3D : public Interpolation3DBase<InternalType, ExternalType> {
 public:
  /* Constructor. Mainly for applying for space.
   * For fixed m, n, l, the class is constructed only once. */
  LinInt3D(int m, int n, int l) {
    // the number of intervals = the number of points - 1
    m--;
    n--;
    l--;
    // initialize *coefficients
    this->coefficients = new InternalType[8 * m * n * l];
    this->m_ = m;
    this->n_ = n;
    this->l_ = l;
  }

  /* The linear form: v = c000 + c100 * x + (c010 + c110 * x) * y +
                         (c001 + c101 * x + (c011 + c111 * x) * y) * z;
   * After calling construct, the 8 coefficients for each grid (m*n*l)
   * are determined, and stored in coefficients[idx],
   * where 0 <= idx < 8*m*n*l. */
  void construct(ExternalType xa[], ExternalType ya[], ExternalType za[],
                 ExternalType ***va) {
    int m = this->m_;
    int n = this->n_;
    int l = this->l_;
    this->xa_ = xa;
    this->ya_ = ya;
    this->za_ = za;
    this->va_ = va;
    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++)
        for (int k = 0; k < l; k++)
          constructGrid(xa, ya, za, va,
                        this->coefficients + 8 * (i * n * l + j * l + k), i, j,
                        k);
  }

  /* The linear form: v = c000 + c100 * x + (c010 + c110 * x) * y +
                         (c001 + c101 * x + (c011 + c111 * x) * y) * z;
   * the numbers in binary form, that is, c100->c[4]
   *       -> cof[0] + cof[4] * x + (cof[2] + cof[6] * x) * y +
   *         (cof[1] + cof[5] * x + (cof[3] + cof[7] * x) * y) * z;
   * First find out the place of (x, y, z) in grid using binary search,
   * then evaluate the polynomial. */
  ExternalType evaluate(ExternalType x, ExternalType y, ExternalType z) {
    int ic = this->binarySearch(this->xa_, this->m_ + 1, x);
    int jc = this->binarySearch(this->ya_, this->n_ + 1, y);
    int kc = this->binarySearch(this->za_, this->l_ + 1, z);
    // printf("%d %d %d \n", ic, jc, kc);
    if (ic > this->m_ - 1) ic = this->m_ - 1;
    if (jc > this->n_ - 1) jc = this->n_ - 1;
    if (kc > this->l_ - 1) kc = this->l_ - 1;
    if (ic < 0 || jc < 0 || kc < 0) {
      if (ic < 0) ic = 0;
      if (jc < 0) jc = 0;
      if (kc < 0) kc = 0;
      return this->evaluateLinear(x, y, z, ic, jc, kc);
    }
    InternalType *cof = this->coefficients +
                        8 * (ic * this->n_ * this->l_ + jc * this->l_ + kc);
    return cof[0] + cof[4] * x + (cof[2] + cof[6] * x) * y +
           (cof[1] + cof[5] * x + (cof[3] + cof[7] * x) * y) * z;
  }

  /* The linear form: v = c000 + c100 * x + (c010 + c110 * x) * y +
                         (c001 + c101 * x + (c011 + c111 * x) * y) * z;
   * the numbers in binary form, that is, c100->c[4]
   *       -> cof[0] + cof[4] * x + (cof[2] + cof[6] * x) * y +
   *         (cof[1] + cof[5] * x + (cof[3] + cof[7] * x) * y) * z;
   * subfunction of construct, compute the coefficients of grid (i, j, k)
   * and store the array in coefficients[8(i*n*l+j*l+k):8(i*n*l+j*l+k)+7].*/
  void constructGrid(ExternalType xa[], ExternalType ya[], ExternalType za[],
                     ExternalType ***va, InternalType *ca, int ib, int jb,
                     int kb) {
    ExternalType x0 = xa[ib], y0 = ya[jb], z0 = za[kb];
    ExternalType dx = xa[ib + 1] - x0, dy = ya[jb + 1] - y0,
                 dz = za[kb + 1] - z0;
    InternalType c000, c001, c010, c011, c100, c101, c110, c111, dt;
    c000 = va[ib][jb][kb];
    c100 = (va[ib + 1][jb][kb] - va[ib][jb][kb]) / dx;
    c010 = (va[ib][jb + 1][kb] - va[ib][jb][kb]) / dy;
    c001 = (va[ib][jb][kb + 1] - va[ib][jb][kb]) / dz;
    c110 = ((va[ib + 1][jb + 1][kb] - va[ib][jb + 1][kb]) / dx - c100) / dy;
    c011 = ((va[ib][jb + 1][kb + 1] - va[ib][jb][kb + 1]) / dy - c010) / dz;
    dt = (va[ib + 1][jb][kb + 1] - va[ib][jb][kb + 1]) / dx;
    c101 = (dt - c100) / dz;
    c111 = (((va[ib + 1][jb + 1][kb + 1] - va[ib][jb + 1][kb + 1]) / dx - dt) /
                dy -
            c110) /
           dz;
    ca[0] = (c000 - c100 * x0) - (c010 - c110 * x0) * y0 -
            ((c001 - c101 * x0) - (c011 - c111 * x0) * y0) * z0;
    ca[4] = (c100 - c110 * y0) - (c101 - c111 * y0) * z0;  // 100
    ca[2] = (c010 - c110 * x0) - (c011 - c111 * x0) * z0;  // 010
    ca[1] = (c001 - c101 * x0) - (c011 - c111 * x0) * y0;  // 001
    ca[6] = c110 - c111 * z0;                              // 110
    ca[5] = c101 - c111 * y0;                              // 101
    ca[3] = c011 - c111 * x0;                              // 011
    ca[7] = c111;
  }

  void printCoefficients() {
    printf("===Debug: Print Coefficients===\n");
    int m = this->m_, n = this->n_, l = this->l_;
    printf("m, n, l: %d %d %d \n", m, n, l);
    InternalType *p = this->coefficients;
    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++)
        for (int k = 0; k < l; k++) {
          printf("Grid (%d, %d, %d) [%.1f, %.1f]x[%.1f, %.1f]x[%.1f, %.1f]: ",
                 i, j, k, this->xa_[i], this->xa_[i + 1], this->ya_[j],
                 this->ya_[j + 1], this->za_[k], this->za_[k + 1]);
          printf(
              "c000 = %.2f, c001 = %.2f, c010 = %.2f, c011 = %.2f, c100 = "
              "%.2f, c101 = %.2f, c110 = %.2f, c111 = %.2f\n",
              p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7]);
          p += 8;
        }
    printf("================================\n");
  }

  ~LinInt3D() { delete[] this->coefficients; }
};

/* QuaInt3D: Quadratic Interpolation for 3 dim
 * */
template <class InternalType = double, class ExternalType = float>
class QuaInt3D : public Interpolation3DBase<InternalType, ExternalType> {
 public:
  /* Constructor. Mainly for applying for space.
   * For fixed m and n and l, the class is constructed only once. */
  QuaInt3D(int m, int n, int l) {
    // the number of intervals = the number of points - 1
    m--;
    n--;
    l--;
    this->m_ = m;
    this->n_ = n;
    this->l_ = l;
    // initialize *coefficients
    this->coefficients = new InternalType[27 * m * n * l];
  }

  /* The quadratic form: z = \sum_{i,j,k} c_{ij} x^i * y^j * z^k.
   * i = 0, 1, 2; j = 0, 1, 2; k = 0, 1, 2.
   * In ternary, i.e. c_121 -> c[16].
   * After calling construct, the 27 coefficients for each grid (m*n*l)
   * are determined, and stored in coefficients[idx],
   * where, 0 <= idx < 27*m*n*l. */
  void construct(ExternalType xa[], ExternalType ya[], ExternalType za[],
                 ExternalType ***va) {
    int m = this->m_;
    int n = this->n_;
    int l = this->l_;
    // initialize *coefficients
    for (int i = 0; i < 27 * m * n * l; i++) this->coefficients[i] = 0.0;

    this->xa_ = xa;
    this->ya_ = ya;
    this->za_ = za;
    this->va_ = va;

    // the inner layer
    for (int i = 0; i < m - 1; i++)
      for (int j = 0; j < n - 1; j++)
        for (int k = 0; k < l - 1; k++)
          constructGrid(xa, ya, za, va,
                        this->coefficients + 27 * (i * n * l + j * l + k), i, j,
                        k);

    // the outer layer
    // i = m-1
    for (int j = 0; j < n; j++)
      for (int k = 0; k < l; k++)
        for (int w = 0; w < 27; w++)
          this->coefficients[27 * ((m - 1) * n * l + j * l + k) + w] =
              this->coefficients[27 * ((m - 2) * n * l + j * l + k) + w];
    // j = n-1
    for (int i = 0; i < m; i++)
      for (int k = 0; k < l; k++)
        for (int w = 0; w < 27; w++)
          this->coefficients[27 * (i * n * l + (n - 1) * l + k) + w] =
              this->coefficients[27 * (i * n * l + (n - 2) * l + k) + w];
    // k = l-1
    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++)
        for (int w = 0; w < 27; w++)
          this->coefficients[27 * (i * n * l + j * l + (l - 1)) + w] =
              this->coefficients[27 * (i * n * l + j * l + (l - 2)) + w];
  }

  /* First find out the place of (x, y, z) in grid using binary search,
   * then evaluate the polynomial. */
  ExternalType evaluate(ExternalType x, ExternalType y, ExternalType z) {
    int ic = this->binarySearch(this->xa_, this->m_ + 1, x);
    int jc = this->binarySearch(this->ya_, this->n_ + 1, y);
    int kc = this->binarySearch(this->za_, this->l_ + 1, z);
    if (ic > this->m_ - 1) ic = this->m_ - 1;
    if (jc > this->n_ - 1) jc = this->n_ - 1;
    if (kc > this->l_ - 1) kc = this->l_ - 1;
    if (ic < 0 || jc < 0 || kc < 0) {
      if (ic < 0) ic = 0;
      if (jc < 0) jc = 0;
      if (kc < 0) kc = 0;
      return this->evaluateLinear(x, y, z, ic, jc, kc);
    }
    // if the point falls in outer layer
    if (ic == this->m_ - 1) ic--;
    if (jc == this->n_ - 1) jc--;
    if (kc == this->l_ - 1) kc--;
    // printf("ic, jc, kc = %d, %d, %d\n", ic, jc, kc);

    InternalType *cof = this->coefficients +
                        27 * (ic * this->n_ * this->l_ + jc * this->l_ + kc);
    ExternalType sum_x = 0, sum_y = 0, sum_z = 0;
    for (int i = 2; i >= 0; i--) {
      sum_y = 0;
      for (int j = 2; j >= 0; j--) {
        sum_z = 0;
        for (int k = 2; k >= 0; k--)
          sum_z = sum_z * (z - this->za_[kc + k]) + cof[i * 9 + j * 3 + k];
        sum_y = sum_y * (y - this->ya_[jc + j]) + sum_z;
      }
      sum_x = sum_x * (x - this->xa_[ic + i]) + sum_y;
    }
    return sum_x;
  }

  /* subfunction of construct, compute the coefficients of grid (i, j, k)
   * and store the array in coefficients[27(i*n*l+j*l+k):27(i*n*l+j*l+k)+26]
   * xa[ib:ib+2], ya[jb:jb+2], za[kb:kb+2] are available data to use. */
  void constructGrid(ExternalType xa[], ExternalType ya[], ExternalType za[],
                     ExternalType ***va, InternalType *ca, int ib, int jb,
                     int kb) {
    // step 1: compute xproduct \prod (x_i-x_j), yproduct, zproduct
    ExternalType xproduct[9], yproduct[9], zproduct[9];
    // xproduct[i][j] = (x_i - x_0)...(x_i - x_{j-1})
    for (int i = 0; i < 3; i++) {
      xproduct[i * 3] = 1;
      for (int j = 1; j < i + 1; j++)
        xproduct[i * 3 + j] =
            xproduct[i * 3 + j - 1] * (xa[ib + i] - xa[ib + j - 1]);
    }
    for (int i = 0; i < 3; i++) {
      yproduct[i * 3] = 1;
      for (int j = 1; j < i + 1; j++)
        yproduct[i * 3 + j] =
            yproduct[i * 3 + j - 1] * (ya[jb + i] - ya[jb + j - 1]);
    }
    for (int i = 0; i < 3; i++) {
      zproduct[i * 3] = 1;
      for (int j = 1; j < i + 1; j++)
        zproduct[i * 3 + j] =
            zproduct[i * 3 + j - 1] * (za[kb + i] - za[kb + j - 1]);
    }

    // step 2: compute a_{pqr} by row
    // temporary vector of \tilde{v}_{ijk}, i = 0, 1, 2, j = 0, 1, 2, k = 0, 1, 2
    InternalType tv[3]; 
    for (int p = 0; p < 3; p++)
      for (int q = 0; q < 3; q++)
        for (int r = 0; r < 3; r++) {
          InternalType res;
          res = va[ib + p][jb + q][kb + r];
          // subtract v_pq0 ... v_pq(r-1), each * zproduct
          for (int k = 0; k < r; k++)
            res -= tv[k] * zproduct[r * 3 + k];
          res /= zproduct[r * 3 + r];
          tv[r] = res;
          // subtract c * xproduct * yproduct
          for (int i = 0; i <= p; i++) {
            InternalType sum_y = 0;
            for (int j = 0; j <= q; j++)
              sum_y += ca[i * 9 + j * 3 + r] * yproduct[q * 3 + j];
            res -= sum_y * xproduct[p * 3 + i];
          }
          res /= xproduct[p * 3 + p] * yproduct[q * 3 + q];
          ca[p * 9 + q * 3 + r] = res;
        }
  }

  void printCoefficients() {
    printf("===Debug: Print Coefficients===\n");
    int m = this->m_, n = this->n_, l = this->l_;
    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++)
        for (int k = 0; k < l; k++) {
          printf(
              "Grid (%d, %d, %d) [%.1f, %.1f] x [%.1f, %.1f] x [%.1f, %.1f]: ",
              i, j, k, this->xa_[i], this->xa_[i + 1], this->ya_[j],
              this->ya_[j + 1], this->za_[k], this->za_[k + 1]);
          for (int k1 = 0; k1 < 3; k1++)
            for (int k2 = 0; k2 < 3; k2++)
              for (int k3 = 0; k3 < 3; k3++)
                printf("a_{%d%d%d} = %.2f; ", k1, k2, k3,
                       this->coefficients[27 * (i * n * l + j * l + k) +
                                          k1 * 9 + k2 * 3 + k3]);
          printf("\n");
        }
    printf("================================\n");
  }

  ~QuaInt3D() { delete[] this->coefficients; }
};

/* SplInt3D: BiCubic Spline Interpolation for 2 dim
 * */
template <class InternalType = double, class ExternalType = float>
class SplInt3D : public Interpolation3DBase<InternalType, ExternalType> {
 private:
  ExternalType **plane_xy;
  InternalType *dy1_, *dy2_, *b_, *T_;

 public:
  /* Constructor. Mainly for applying for space.
   * For fixed m, n and l, the class is constructed only once. */
  SplInt3D(int m, int n, int l) {
    this->m_ = m - 1;
    this->n_ = n - 1;
    this->l_ = l - 1;
    // initialize *coefficients
    this->coefficients = new InternalType[m * n * l];
    // in evaluate
    plane_xy = new ExternalType *[m];
    for (int i = 0; i < m; i++) plane_xy[i] = new ExternalType[n];
    // in construct 1D
    // first-order & second-order difference
    int N = l;
    dy1_ = new InternalType[N];
    dy2_ = new InternalType[N];
    b_ = new InternalType[N];
    T_ = new InternalType[N];
  }

  /* After calling construct, the second-order derivatives M_i
   * (0 <= i < l)
   * for each pair (x, y) are determined, and stored in
   * coefficients[(i * n + j) * l : (i * n + j) * l + l - 1] */
  void construct(ExternalType xa[], ExternalType ya[], ExternalType za[],
                 ExternalType ***va) {
    int m = this->m_ + 1;
    int n = this->n_ + 1;
    int l = this->l_ + 1;
    this->xa_ = xa;
    this->ya_ = ya;
    this->za_ = za;
    this->va_ = va;
    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++)
        construct1D(za, va[i][j], this->coefficients + (i * n + j) * l, l);
  }

  /* First find out the place of (x, y, z) in grid using binary search,
   * then evaluate the polynomial. */
  ExternalType evaluate(ExternalType x, ExternalType y, ExternalType z) {
    int ib = this->binarySearch(this->xa_, this->m_ + 1, x);
    int jb = this->binarySearch(this->ya_, this->n_ + 1, y);
    int kb = this->binarySearch(this->za_, this->l_ + 1, z);
    if (ib < 0 || jb < 0 || kb < 0 || ib > this->m_ - 1 || jb > this->n_ - 1 ||
        kb > this->l_ - 1) {
      if (ib < 0) ib = 0;
      if (jb < 0) jb = 0;
      if (kb < 0) kb = 0;
      if (ib > this->m_ - 1) ib = this->m_ - 1;
      if (jb > this->n_ - 1) jb = this->n_ - 1;
      if (kb > this->l_ - 1) kb = this->l_ - 1;
      return this->evaluateLinear(x, y, z, ib, jb, kb);
    }
    int m = this->m_ + 1;
    int n = this->n_ + 1;
    int l = this->l_ + 1;
    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++)
        plane_xy[i][j] =
            evaluate1D(this->za_, this->va_[i][j],
                       this->coefficients + (i * n + j) * l, l, z, kb);
    // Given plane, construct 2D
    SplInt2D<InternalType, ExternalType> s =
        SplInt2D<InternalType, ExternalType>(m, n);
    s.construct(this->xa_, this->ya_, plane_xy);
    return s.evaluate(x, y);
  }

  /* Construct 1D spline, the arrays are of size n*/
  void construct1D(ExternalType ya[], ExternalType za[], InternalType *ca,
                   int n) {
    // first-order & second-order difference
    InternalType *dy1 = dy1_;
    InternalType *dy2 = dy2_;
    for (int i = 0; i < n - 1; i++) dy1[i] = ya[i + 1] - ya[i];
    for (int i = 0; i < n - 2; i++) dy2[i] = ya[i + 2] - ya[i];
    InternalType *b = b_;
    for (int i = 0; i < n - 1; i++) b[i] = (za[i + 1] - za[i]) / dy1[i];
    for (int i = 0; i < n - 2; i++) b[i] = 6 * (b[i + 1] - b[i]) / dy2[i];

    // solve Tx = b. T stores diagonal of T
    // TODO: Use LAPACK library to solve triangular systems?
    InternalType *T = T_;
    T[0] = 2;
    for (int i = 1; i < n - 2; i++) {
      InternalType scale = dy1[i] / dy2[i] / T[i - 1];
      T[i] = 2 - dy1[i] / dy2[i - 1] * scale;
      b[i] = b[i] - b[i - 1] * scale;
    }
    ca[0] = 0;
    ca[n - 1] = 0;
    if (n > 2) ca[n - 2] = b[n - 3] / T[n - 3];
    for (int i = n - 3; i > 0; i--)
      ca[i] = (b[i - 1] - dy1[i] / dy2[i - 1] * ca[i + 1]) / T[i - 1];
  }

  /* f = c0 + c1*x + c2*x*x + c3*x*x*x
   *   = c0 + x*(c1 + x*(c2 + x*c3))
   * */
  InternalType evaluate1D(ExternalType xa[], ExternalType za[],
                          InternalType *ca, int n, ExternalType x, int i) {
    x -= xa[i];
    InternalType c0, c1, c2, c3, dx = xa[i + 1] - xa[i];
    c0 = za[i];
    c1 = (za[i + 1] - za[i]) / dx - (ca[i + 1] / 6 + ca[i] / 3) * dx;
    c2 = ca[i] / 2;
    c3 = (ca[i + 1] - ca[i]) / dx / 6;
    return c0 + x * (c1 + x * (c2 + x * c3));
  }

  void printCoefficients() {
    printf("===Debug: Print Coefficients===\n");
    int m = this->m_ + 1, n = this->n_ + 1, l = this->l_ + 1;
    for (int i = 0; i < m; i++)
      for (int j = 0; j < n; j++)
        for (int k = 0; k < l; k++)
          printf(
              "(%d, %d, %d): x = %.3f, y = %.3f, z = %.3f, v = %.3f, ca = %.3f\n",
              i, j, k, this->xa_[i], this->ya_[j], this->za_[k],
              this->va_[i][j][k], this->coefficients[i * n * l + j * l + k]);
    printf("================================\n");
  }

  ~SplInt3D() {
    delete[] this->coefficients;
    for (int i = 0; i < this->m_ + 1; i++) delete[] plane_xy[i];
    delete[] plane_xy;
    delete[] dy1_, dy2_, b_, T_;
  }
};

#endif